function mmrpz = fun_mrpz(Par,agg,mz,ml)
% Returns the marginal revenue product of productivity (MRPZ) with given 
% productivity and labor allocation
% Output:
%   mmrpz - matrix of marginal revenue product of productivity with 4 rows
%   and column number consistent with those of mz and ml
% Input:
%   Par - parameter structure
%   agg - 1 by 3 vector of equilibrium values of aggregate variables
%   mz - productivity matrix, 2 rows, same num of columns as ml
%   ml - matrix of labor allocation, 4 rows, same num of columns as mz
    sigma = Par.sigma;
    rho = Par.rho;
    tau = Par.tau;
    moutput = kron(mz,[1;1]).*ml./[1;tau;tau;1];
    sector_output = [(moutput(1,:).^((rho-1)/rho)+moutput(3,:).^((rho-1)/rho)).^(rho/(rho-1));(moutput(2,:).^((rho-1)/rho)+moutput(4,:).^((rho-1)/rho)).^(rho/(rho-1))];
    sector_price = ([agg(1);agg(2)]./sector_output).^(1/sigma);
    mprice = repmat(sector_price.*sector_output.^(1/rho),2,1)./moutput.^(1/rho)./[1;tau;tau;1];
    momega = (moutput./repmat(sector_output,2,1)).^((rho-1)/rho);
    mepsilon = 1./(momega./sigma+(1-momega)./rho);
    mmarkup = mepsilon./(mepsilon-1);
    mmrpz = ml.*mprice./mmarkup;
end